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Progress in the understanding of quantum critical properties of itinerant electrons has been hin¬ 
dered by the lack of effective models which are amenable to controlled analytical and numerically 
exact calculations. Here we establish that the disorder driven semimetal to metal quantum phase 
transition of three dimensional massless Dirac fermions could serve as a paradigmatic toy model for 
studying itinerant quantum criticality, which is solved in this work by exact numerical and approx¬ 
imate field theoretic calculations. As a result, we establish the robust existence of a non-Gaussian 
universality class, and also construct the relevant low energy effective field theory that could guide 
the understanding of quantum critical scaling for many strange metals. Using the kernel polynomial 
method (KPM), we provide numerical results for the calculated dynamical exponent (z) and corre¬ 
lation length exponent (u) for the disorder-driven semimetal (SM) to diffusive metal (DM) quantum 
phase transition at the Dirac point for several types of disorder, establishing its universal nature 
and obtaining the numerical scaling functions in agreement with our field theoretical analysis. 
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I. INTRODUCTION 

Phase transitions are ubiquitous in the natural world 
ranging from the solidification of water to the thermal 
suppression of magnetism in iron or superconductivity in 
metals. Through the construction of simplified effective 
“toy” models, such as the Ising model to describe clas¬ 
sical magnets [I|, generic theories can be developed to 
gain invaluable physical insights and explain the univer¬ 
sal features of experimental data across vastly different 
systems. Due to the quantum mechanical zero-point mo¬ 
tion, phase transitions can also occur at absolute zero 
temperature 0 driven entirely by quantum (rather than 
thermal) fluctuations, which can be accessed by tuning a 
non-thermal control parameter g, such as disorder, dop¬ 
ing, magnetic held, or pressure. When such a quantum 
phase transition (QPT) is continuous, the quantum crit¬ 
ical point (QCP) (located at a critical coupling gc akin 
to the critical temperature for thermal phase transitions) 
exhibits scale invariance, and the universal scaling prop¬ 
erties of physical quantities are governed by the divergent 
length (^i) and time (^t) scales 

(i-\9-9cr, 

where the critical exponents v and z characterize the 
correlations in space and time respectively. These no¬ 
tions have been put on a solid theoretical foundation, 
e.g. by addiM non-trivial quantum dynamics into the 
Ising model Q, which has allowed for a simplihed set¬ 
ting to understand magnetic quantum phase transitions 
in general. The critical exponents v and z dehne the 
universality class of the quantum phase transition, and 
critical scaling functions connecting various physical pa¬ 
rameters define the nature of the QCP. 

Precisely at the QCP {g = gc) and at T = 0, one 
has true scale invariance as the characteristic correlation 


lengths and are infinite. Even though the QCP 
occurs strictly at zero temperature, its effects on physical 
properties manifest over a broad range of temperatures 
(T) and energies {E), such that 

E, ksT > 

where ks is the Boltzmann constant. This critical re¬ 
gion, known as the quantum critical fan (see Fig. 1), 
is naturally accessible in the laboratory and is widely 
studied [3 t 9- The generic magnetic QCPs for most in¬ 
sulating systems can be simply understood in terms of 
the Landau-Ginzburg-Wilson order parameter theory in 
(d-|-z) dimensions, since there is only one gapless bosonic 
degree of freedom, namely the order parameter. By con¬ 
trast, at the the metallic or itinerant QCPs there are 
at least two sets of gapless degree of freedom, (i) the 
itinerant fermions and (ii) the bosonic order parameter, 
whose mass has been tuned to zero (at g = gc). Thus, to 
construct a successful description of a metallic or itiner¬ 
ant QCP, it is necessary to address both types of gapless 
excitations on an equal footing, which is generally a chal¬ 
lenging task. 

In this work, we show that the disorder driven 
semimetal to diffusive metal (SM-DM) QCP of three- 
dimensional Dirac fermions is a quintessential toy model 
for gaining qualitatively new insights into the general 
nature of itinerant quantum criticality. For a random 
scalayjotential, this problem has been studied by mean- 
field 0 and perturbative field theoretical calculations [1- 
[l^ . These analytical calculations show that above the 
lower critical dimensionality d = 2, the average density 
of states (DOS) at zero energy vanishes inside the SM for 
sufficiently weak disorder, and it becomes finite beyond 
a critical disorder strength leading to a diffusive metal 
(DM). There is no such finite disorder phase transition 
in two dimensions (e.g. graphene), where the system be- 
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FIG. 1: (Color Online) Numerically calculated finite temper¬ 
ature (T) phase diagram as a function of disorder {W) com¬ 
puted for a linear system size L = 110 (with t, the band hop¬ 
ping energy, setting the energy scale). The Dirac semimetal 
(SM), the diffusive metal (DM) and the QCP separating them 
only exist at zero temperature and there is no finite temper¬ 
ature phase transition. The squares mark the cross over out 
of the SM phase (green) into the quantum critical (QC) re¬ 
gion (orange), which is anchored by the QCP at the critical 
disorder value Wc/t = 2.55 ± 0.05 (star). This cross over is 
determined by where the specific heat fails to be described 
by Cv ^ ■ The circles mark the cross over out of the DM 

phase (blue) into the QC region, which is marked by where the 
specific heat fails to be described by Cv ~ T. The dashed 
lines connecting the squares/circles are a guide to the eye. 
The dashed line at A/f marks the high energy cut off above 
which the low energy continuum description in terms of mass¬ 
less Dirac fermions is inapplicable. In a narrow region inside 
the quantum critical regime the specific heat manifests non- 
Fermi-liquid behavior Cv ~ which gradually crosses over 
to the power law behaviors of SM and DM. The dashed hor¬ 
izontal line Ep/t is a schematic for the (generally unknown) 
Fermi energy associated with any residual doping that drives 
the system slightly away from the precise Dirac point, indicat¬ 
ing the low-energy cut off above which the zero-doping Dirac 
point theory is applicable. Note that as long as Ep is not 
too large, quantum criticality manifests itself in the orange 
region of the phase diagram, even though the QCP itself (the 
star at Wc and T = 0) is ‘hidden’ by the residual doping. At 
temperature and energy scales smaller than Ep, the system 
behaves as a conventional metal or diffusive Fermi liquid. In 
a nominally undoped system, the inequality Ep applies, 
and quantum criticality is observable for kpT > Ep. 

comes a diffusive metal already for infinitesimal disorder. 
Through a, d = 2 -\- e expansion of the critical properties, 
it is found that the leading order results for the critical 
exponents in d = 3 are z = 3/2 and iz = I, for both po¬ 
tential and axial disorder, while there should be no such 
transition for mass disorder Q. Recently, numerical cal¬ 
culations have also come to bear on the problem through 
the study of a three dimensional disordered topological 
insulato r fl^ . [3, a three-dimensional layered Chern in¬ 
sulator [2, a Weyl semimetal lla-fla . and the phase 
diagram of Dirac and Weyl semimetals. For 
the case of a single Weyl cone which can be realized on 


the surface of a four dimensional topological insulator, 
the critical exponents have been numerically obtained to 
high accuracy M- For avoiding any misunderstanding, 
we mention that there is a second quantum phase transi¬ 
tion for stronger disorder {Wi » Wc than considered 
in Fig. [IJ, where the diffusive metal undergoes an An¬ 
derson localization. In this work we do not address the 
Anderson localization transition at all and exclusively fo¬ 
cus on the SM-DM QPT restricting ourselves to disorder 
strengths below the threshold for the Anderson transi¬ 
tion (1^ . 

This itinerant QPT is addressed in our work through 
numerically exact calculations on sufficiently large sys¬ 
tem sizes (~ 10® lattice points) and approximate field 
theoretic analysis, which we use to establish: 

1. The numerical value of z is universal, being in¬ 
dependent of the disorder type and the details of 
their probability distribution (box versus a Gaus¬ 
sian distribution) and we obtain numerically z = 
1.46 ± 0.05 (in good agreement with other numeri¬ 
cal studies on different models 3E3). For the 
correlation length exponent we cannot pinpoint a 
universal value and we find a value of v that lies in 
the range ^ 0.9 — 1.5. (Our numerical technique, 
while providing a fairly precise value of z, is rela¬ 
tively imprecise in obtaining iz because of the uncer¬ 
tainties in the precise determination of the critical 
disorder strength Wc which turns out to be crucial 
in calculating u, but not z.) 

2. The existence of a broad quantum critical fan at 
finite temperatures (see Fig. 1), where the fan oc¬ 
cupies a large parameter space. 

3. We find considerable agreement between the nu¬ 
merically exact crossover scaling functions and 
their approximate analytical forms (Appendix A) 
for the average density of states and the specific 
heat. The approximate analytical calculations are 
based on one loop renormalization group analysis 
which predicts z = 3/2 and v = 1. 

4. We also construct an order parameter field the¬ 
ory coupled to itinerant fermions, which illustrates 
how the notion of well defined quasiparticle exci¬ 
tations can break down in the quantum critical 
regime. The microscopic model we solve is non¬ 
interacting (with disorder as the tuning parameter 
for the QCP), and the sample to sample fluctu¬ 
ations of disorder give rise to effective electronic 
interactions of a statistical nature. Even though 
this does not pertain to correlated electronic sys¬ 
tems per se, the theoretical concepts, regarding the 
itinerant quantum critical scaling properties that 
we are able to establish, should be generally appli¬ 
cable to metallic quantum critical points. In par¬ 
ticular, our theoretical finding and the numerical 
verification of the non-Fermi-liquid behavior and 
crossover scaling functions in the quantum critical 
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regime should have generic qualitative applicabil¬ 
ity to itinerant quantum critical systems. There¬ 
fore, this model may effectively serve as an “Ising 
model” for itinerant quantum criticality for future 
theoretical work. 

5. Our theoretical results are directly pertinent for de¬ 
scribing the recently discovered gapless Dirac semi¬ 
conductors or semimetals | 2 ll - l^ . where the va¬ 
lence and conduction bands touch linearly at iso¬ 
lated points in the Brillioun zone. In particular, we 
expect that the quantum critical fan established in 
this work will be experimentally observable in Dirac 
semimetals with a very low carrier concentration. 
We therefore pro pose NaaBi, which has a small 
Fermi energy (24| . as one of the most promising 
candidate materials for realizing the predicted crit¬ 
ical properties. Our results can also be relevant for 
the three dimensional topological quantum phase 
transitions between a topological and band insula¬ 
tor in the presence of disorder (26l - [30l| . Although 
our quantum critical results strictly apply only for 
the intrinsic (i.e. undoped) Dirac semimetal, the 
results should remain valid even for finite doping 
as long as the doping density is low enough so that 
the associated chemical potential is smaller than 
the temperature [3l|. 

In our numerically calculated phase diagram 
(Fig.[T|) we depict the SM, DM, and the critical fan 
regions with green, blue, and orange colors respec¬ 
tively whereas the (bottom) dashed horizontal line 
schematically indicates the Fermi level (or chemi¬ 
cal potential) associated with any residual (unin¬ 
tentional and hence unknown) doping in the sys¬ 
tem providing the low-energy cut off for the quan¬ 
tum critical theory. The high-energy cut off is in¬ 
dicated by the top horizontal line which is associ¬ 
ated with the usual ultraviolet cut off for studying 
Dirac physics on a lattice. We note that the invari¬ 
able presence of some residual doping in the sys¬ 
tem producing the low-energy cut-off Ep in Fig. [1] 
acts to ‘hide’ the QCP by introducing an incipient 
(unintentional) metallic phase. Although the QCP 
itself is inaccessible and hidden, its effect persists 
in the quantum critical regime (the orange region 
above the Ep line in Fig.[T|) provided that the dop¬ 
ing is not too high. This situation is not dissimilar 
to many examples of itinerant quantum criticality 
in correlated metallic systems. 

The remainder of the paper is organized as follows: In 
section El we introduce the lattice model and its contin¬ 
uum limit. In section uni we present detailed numerical 
results, which is followed by the construction of the or¬ 
der parameter field theory in section HVl and we conclude 
in section [V] In Appendix we theoretically determine 
the general scaling formulas used to analyze the numer¬ 
ical data. Finally additional numerical details for our 
calculations are provided in Appendix |B] 


II. LATTICE MODEL AND CONTINUUM 
LIMIT 


We begin by introducing the lattice model and its con¬ 
tinuum limit that will be the focus of this work. We 
consider non-interacting electrons hopping on a simple 
cubic lattice with periodic boundary conditions in the 
presence of a random on site potential. The Hamiltonian 
is 

H “jV’r+e, + H.c) + ^ P {y)iI}IAwiPv (1) 

rj r 


where is a four component spinor. For non¬ 
superconducting systems with conserved electric charge 
we can choose ip'^ = (cr,+,'|-, Cr,-,-!-, Cr,+, 4 ,, Cr,-, 4 ,) which is 
composed of an electron annihilation operator Cr,r,s at 
site r, with parity r = ±, and spin-projections s =t / i- 
We work in the Dirac representation and therefore the 
matrices are 


a,- = 


0 cr, 

CTj 0 



/3 


1 0 

0 -1 


( 2 ) 


where aj denotes the Pauli matrices and 1 denotes the 
2x2 identity matrix. Each site is labeled by the vector 
r and the nearest neighbors are connected by the unit 
vectors ej, with j = 1, 2,3. The strength of hopping be¬ 
tween neighboring sites is determined by t, and the ran¬ 
dom potential is P(r). (We often use t = 1 in the rest of 
this paper using the hopping as the unit of energy in the 
problem.) To study the universality of the transition we 
consider two types of random distributions for V (r) : (i) 
a box distribution (BD), i.e. a randomly distributed vari¬ 
able between [—IF/2, VF/2], (ii) a Gaussian distribution 
(GD) with zero mean and variance W^. The type of dis¬ 
order is specified by the Aw matrix, we consider three 
different types of disorder: ( 1 ) axial disorder A 5 = 75 , 
( 2 ) mass disorder Am = /3, and (3) potential disorder 
Ap = Iix 4 (where denotes the four by four identity 
matrix). Physically, each type of disorder can naturally 
occur in experimental systems without any control over 
which one may appear or produce the dominant effect, 
and it is not unreasonable to assume that the QPT tuned 
by each type of disorder belongs to a different universal¬ 
ity class (which turns out not to be the case here as we 
will see). 

In addition to the global U(l) symmetry for total 
number conservation, the tight binding model also has 
a continuous U(l) chiral symmetry, as the Hamiltonian 
commutes with ''/’rTs'/’r, where 75 = iQ:ia 2 Q; 3 . This 
Hamiltonian serves as a model for a topological Dirac 
semimetal with eight Dirac cones at the high symme¬ 
try points of the cubic Brillouin zone, and the excita¬ 
tions around different cones may be coupled by inter¬ 
valley scattering due to disorder. As a result of the 
U(l) chiral symmetry we expect the axial and poten¬ 
tial disorders to have identical critical properties, man¬ 
ifesting universality. An important question is whether 
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the universality of the SM-DM QCP remains unaffected 
when the C/(l) chiral symmetry is absent. We address 
this question by studying the effects of mass disorder, 
which reduces the U(l) chiral symmetry down to a Zi 
chiral or particle-hole symmetry, described by the rela¬ 
tion {i?, ''/^r*/375'0r} = 0. While potential disorder is 

only relevant for non-superconducting systems, axial and 
mass disorder can naturally appear for superconducting 
systems. Due to the discrete particle-hole symmetry, the 
mass disorder problem can appear for class AIII or Dill 
(depending on the representation of the spinor). By con¬ 
trast, the axial disorder displays the discrete particle hole 
symmetry with respect to both and any arbitrary 

linear combination of them. Thus, the axial disorder 
model is shared by both AIII and Dill Altland-Zirnbauer 
classes. One of our completely unanticipated interesting 
new findings is that the average density of states at zero 
energy (p(0)) across the SM-DM QCP is independent of 
the type of disorder (potential, mass, or axial) driving 
the transition, see Fig. [5] (a). As p(0) acts as the order 
parameter of the SM-DM transition, this implies same 
power law dependence of p(0) for potential, axial, and 
mass disorder. Consequently, the SM-DM QPTs driven 
by potential, mass, and axial disorders belong to the same 
universality class. We establish this universality of the 
disorder-tuned SM-DM transitions through nonperturba- 
tive exact numerical calculations. 

In the absence of disorder, the model reduces to the 
well known Dirac Hamiltonian on a lattice, with a dis¬ 
persion 


Eo{k) = ±t^ /^sin(fcj) 


( 3 ) 


Expanding the dispersion in the low energy limit gives 
rise to eight Dirac cones with linearly dispersing exci¬ 
tations around eight high symmetry points of the cubic 
Brillouin zone. Considering only the long wavelength de¬ 
grees of freedom (which are the most dominant at a sec¬ 
ond order phase transition) in the presence of disorder 
we arrive at the following continuum Hamiltonian as the 
effective theory 

'^D = ipi {-ivajdj Zilxa + V{x)Aw 0 Bxa) i’a, (4) 

where Ixa is an eight by eight identity matrix in the val¬ 
ley space, Bxy is the intervalley mixing matrix, and now 
the Dirac spinor carries a valley or flavor index A. In 
the absence of intervalley scattering there is an SU{8) 
flavor symmetry. The intervalley scattering due to disor¬ 
der breaks this flavor symmetry. We mention, however, 
that for long-ranged disorder, e.g.. Coulomb impurities, 
such a disorder-induced valley symmetry breaking may 
be weak and thus operational only at very low temper¬ 
atures, which appears to be the situation in (the two- 
dimensional Dirac material) graphene. In such a situa¬ 
tion, the independent Dirac cone approximation that ig¬ 
nores intervalley scattering effects is valid down to very 
low temperatures. 


For simplicity of the analysis, the analytic calcula¬ 
tions are performed with the independent Dirac cone ap¬ 
proximation, while neglecting intervalley scattering (i.e. 
Bxa —>■ Ixa)- We choose a Gaussian white noise disorder 
distribution with a variance A, and define the dimension¬ 
less disorder strength A = AA'^“^Hd/[(27r)^u^] where A 
is the large momentum cut off and = 27r‘^/^/r(d/2) 
is the area of the unit sphere in d dimensions, and F(a;) 
is the gamma function. Note that for analytic calcula¬ 
tions we refer to the reduced distance from the QCP as 
i5 = |A — Ac|/Ac, while for the numerical work we de¬ 
note the reduced distance as 5 = \W — Wc\/Wc- To avoid 
confusion we explicitly mention our definitions when we 
use 5. The one loop calculation (or to the lowest order 
in e) for potential or axial disorder yields the same beta 
function and the SM-DM QPT is controlled by the 
QCP located at Ac = {d — 2)/2 = e/2 with z = 1 + e/2 
and IX = l/{d — 2). For d = 3, this leads to z/ = I, and 
z = 3/2. The implications of one loop calculations for 
determining quantum critical scaling properties are dis¬ 
cussed in detail in Appendix]^ By contrast, the one-loop 
beta function shows intravalley mass disorder to be an 
irrelevant perturbation and does not predict any QPT. 
However, through our nonperturbative numerical calcu¬ 
lations we show that the mass disorder with intervalley 
scattering actually drives a SM-DM QPT (see Fig. [5] (a)), 
in contrast to the one loop perturbative result with only 
intravalley scattering. 

The important question is whether the results derived 
at the one loop level for potential and axial disorders in 
the single-cone approximation remains valid in the pres¬ 
ence of intervalley scattering, i.e., how robust the analyti¬ 
cally obtained values of the critical exponents and scaling 
functions in the general situation are where the precise 
theoretical assumptions necessary for obtaining the ana¬ 
lytical results may no longer apply. To answer these ques¬ 
tions non-perturbatively, we have performed detailed nu¬ 
merical calculations on the tight binding model in Eq. (HD 
using large system sizes, which we present in the next 
section. 


III. NUMERICAL RESULTS 

We are primarily interested in obtaining the disorder 
averaged DOS. As shown below the DOS at zero energy 
serves as the order parameter for the SM-DM QPT. The 
DOS is defined as 


p{E,L) = ^^SiE-E,), (5) 

i=l 

where L is the linear size of the system with a vol¬ 
ume V = D = AV is the total number of states 
of the model, and Ei corresponds to each eigenstate. 
We calculate the DOS using the numeric ally exact ker¬ 
nel polynomial method (KPM) (see Ref. for explicit 
details), which allows us to reach sufficiently large sys- 
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FIG. 2: (Color Online) Numerically (KPM)-calculated DOS as a function of disorder strength for a linear system size L = 110. 
(a) DOS at zero energy p(0) as a function of disorder strength for each type of disorder matrix Aw, which shows that axial, 
mass, and potential disorder yield identical results. (Inset) Corresponding calculation for axial disorder in two dimensions [i.e. 
setting the hopping in the ^-direction to zero in Eq. 0] clearly establishing that there is no stable SM phase and p(0) is finite 
even for an infinitesimal amount of disorder, (b) DOS for axial disorder as a function of energy passing through the QCP with 
the arrow indicating increasing disorder strength. Specific heat (solid lines) with axial disorder for W < Wc (c) and W > Wc 
(d), the fits (dashed lines) are to the leading Cv ~ in the SM and Cv ~ T in the DM. The temperatures at which the fits 
no longer apply determine the cross over boundary shown in Fig. [T] 


tem sizes (also see Appendix IbI) . We calculate the spe¬ 
cific heat Cv{T) = d{E)/dT by numerically integrating 
the density of states using the following formula for non¬ 
interacting fermions 


Cv 



dE 


p{e)e^ 

4cosh(/3£’/2)2 ’ 


( 6 ) 


where /3 = 1/T is the inverse temperature and we are 
working in units where fcs = 1. 


A. Method 

Central to KPM [32|, we expand the density of states 
in terms of Chebyshev polynomials Tn{x) and truncate 
the expansion at order n = N^,. (The dependence of 
the numerical convergence on the parameter N^. must 
be checked explicitly- see Appendix [B] for details.) 


As Chebyshev polynomials are only dehned on an 
interval [—1,1], we have to rescale the Hamiltonian to 
ensure that its eigenvalues fall within the corresponding 
range. This can be done from a simple transformation 
E[ = aEl' + b where a and b are related to the maximum 
and minimum eigenvalues, which we estimate using the 
Lanczos method. In short, the KPM reduces the prob¬ 
lem of diagonalizing the Hamiltonian into calculating 
the coefficients of the expansion = Tr[r„(iT')], which 
can be done using only matrix-vector multiplication, 
and the trace can be evaluated stochastically. For the 
sparse Hamiltonian matrix considered here, this can be 
done very efficiently, and therefore, KPM is capable of 
handling system sizes much larger than what can be 
done by direct diagonalization for our purpose. As is well 
known from the analysis of Fourier series, truncating a 
series expansion can lead to artificial oscillations (Gibb’s 
oscillations) in the calculation. We hlter out these 
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FIG. 3: (Color Online) Determination of the critical exponent 2 for potential disorder with the box distribution and L = 130. (a) 
DOS as a function of energy, and (b) the specific heat as a function of temperature, both obtained at the QCP, Wc = 2.55±0.05. 


spurious effects by using the Jackson Kernel [l^. In 
Appendix |B] we discuss in detail the convergence of the 
KPM and how to choose appropriate values for Nc- For 
the results presented in Figs.[T]and[ll we have calculated 
Nc = 1028 moments such that the DOS is smooth, while 
considering a lattice size L = 110 (i.e. a system size 
V = 110^, which would be unthinkable from the exact 
diagonalization perspective). As the density of states 
is self averaging (which we have explicitly checked), 
we only have to average over different realizations of 
the disorder to reduce the noise in the calculation and 
obtain a smooth DOS. For the critical exponents, we 
consider Nc = 4096 and focus on a two component 
model (i.e. replacing all the aj matrices in Eq. (ED 
with Pauli matrices aj) after isolating one of the two 
degenerate eigenvalues due to the axial symmetry. In 
this case, we consider system sizes ranging from L = 60 
up to L = 130, and in order to completely eliminate 
any statistical fluctuations we perform 100 disorder 
averages. The scaling formulas used for the analysis of 
the numerical data are described in Appendix 


B. Potential, Axial, and Mass Disorder 

As shown in Fig. [Hb), the DOS inside the SM phase 
indeed scales as p(-E) oc |i?p, making the SM an in¬ 
compressible, gapless state, which remains stable up to 
a critical strength of disorder. For all three types of dis¬ 
order (potential, axial, and mass), p(0) becomes finite 
after passing through a QCP at a finite disorder strength 
Wc ~ 2.55t, the putative SM-DM transition driven by 
increasing disorder. Remarkably, the microscopic criti¬ 
cal coupling for all three types of disorder are identical 
within our numerical accuracy, as clearly demonstrated 
in Fig. 2(a). This is a rather amazing unanticipated 
result since universality in critical phenomena usually 


refers to the universality of the critical exponents (and 
the scaling functions in dimensionless units), but not to 
the critical coupling strength (or the critical tempera¬ 
ture Tc in thermodynamic phase transitions). The fact 
that the mass disorder leads to the QPT highlights the 
importance of intervalley scattering. At the same time 
the equal strength of microscopic critical couplings for all 
three types of disorder can not be addressed within any 
effective low energy theory and we have no theoretical 
explanation for this finding, since the critical coupling 
strength is not known to be a universal quantity, in con¬ 
trast to critical exponents and scaling functions, in any 
quantum (or classical) critical phenomenon. More work 
would be necessary, which is well beyond the scope of the 
current work, to understand this unexpected ‘universality 
in critical coupling strength’ which we have numerically 
discovered in this problem. 

By contrast, the DOS for a two dimensional Dirac 
SM becomes finite for an infinitesimally weak disorder 
strength as shown in the inset of Fig. 2(a), which agrees 
with the field theoretic predictions that d = 2 is the 
lower critical dimension for the disorder-tuned SM-DM 
QPT. The behavior of the three dimensional SM is also 
clearly distinct from the effects of disorder in compress¬ 
ible normal metals, which for an infinitesimal amount 
of disorder converts into a DM (from a ballistic metal 
with perfect conductivity). We have explicitly shown 
in Ref. [l^ that the DM phase (in the present model) 
is not Anderson localized, and only for a large disorder 
strength Wi/t k. 8.8 (for the box distribution) does the 
DM undergo Anderson localization, similar to what hap¬ 
pens in three dimensional conventional disordered met¬ 
als [33 - 1^ . Thus, within the present calculation we find 
the SM (for W < Wc) and the DM (for W > Wc) phases 
to be stable over finite ranges of disorder in the undoped 
three dimensional Dirac system and we use the specific 
heat to determine the cross over energy scales for each 
phase as shown in Figs. [5] (c) and (d). 







(b) 



FIG. 4: (Color Online) Determination of the critical exponent v for potential disorder with the box distribution [(a) and (c)] 
and the Gaussian distribution [(b) and (d)]. (a) and (b) show the numerical DOS as a function of the distance to the QCP. 
The dashed lines are power law fits to extract v, the thick (thin) dashed lines are fits closer to (further from) the QCP. (c) and 
(d) show the numerical finite size scaling and data collapse of p(0, L), the labels for L are shared between (c) and (d). (Insets) 
The local-linearity function [^, which yields the best scaling collapse at its minimum value. Here 5 — \W — 1141/114 is the 
usual dimensionless tuning parameter defining the QPT. 


C. Universal scaling and critical exponents 


From our numerical results we now show that the den¬ 
sity of states and specific heat obey universal scaling 
forms in the vicinity of the QCP, which enables us to de¬ 
termine the critical exponents characterizing the univer¬ 
sality class (see Appendix]^ for the details regarding the 
scaling ansatz). Assuming that the QCP is non-Gaussian 
(i.e. non-mean field), we apply the scaling hypothesis to 
the density of states and specific heat implying that they 
must obey the scaling forms 


p{E,L) = ,5(^-")"7^(|^;|(5-"^Ll/"(5) (7) 

Cv{T) = (8) 


where we have introduced the reduced distance to the 
QCP, 6 = \W — Wc|/bFc, and TZ and C are two unknown 


scaling functions, which are related through 

'Y'£ — L/Z P + OO 

C{TS-''\L^^^S) = —— / duu^ cosh-2(u/2) 

^ J —oo 


xn{\u\T5-''\L^/''5). 


(9) 


Based on the numerical calculations presented in Fig. [2] 
(a), we have established that the QCP driven by axial 
or mass or potential disorder between the SM and the 
DM falls within the same universality class (within our 
numerical accuracy). Therefore, we conclude that irre¬ 
spective of the underlying chiral symmetry being present 
or absent, the disorder driven SM-DM QPT of three- 
dimensional massless Dirac fermions exhibits a manifest 
universality, which is one of our main results. 

We determine the location of the critical point by ex¬ 
trapolating p(0) to zero from the DM phase, which yields 
Wc/t = 2.55 ±0.05 for box disorder (and for all three dif¬ 
ferent types of disorder) and Wc/t = 0.60±0.03 for Cans- 
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(b) 






FIG. 5: (Color Online) Scaling collapse of the DOS and the specific heat in the quantum critical fan 5'^^ <C Elt,Tlt <C A/t 
(with (5 = \W — Wc|/Wc) establishing that both exhibit single parameter scaling as a function of the distance to the QCP given 
by Eqs. 0 and (121). The dashed line is the approximate crossover functions analytically determined from the one loop beta 
functions (see Appendix |3; (a) and (c) are for W < Wc, and (b) and (d) are for W > Wc- We stress that the dashed lines 
are not a fit, and only the non-universal amplitudes are adjusted. Despite the significant difference between the numerically 
obtained value for v and the one loop result {v — 1), we find considerable agreement between the analytical and numerical 
crossover functions. We are using 2 : = 1.46 and the value of v from the finite size scaling for the collapse of the numerical 
data vl ~ 1.46. These results are for the axial disorder using the box distribution, but the results for other distributions and 
disorder are very similar. 


sian potential disorder. Thus, rather curiously, while our 
numerically obtained critical disorder Wc is the same for 
the three types of disorder we consider (i.e. potential, 
mass, axial), it is not for the different forms of the disor¬ 
der distribution (box or Gaussian) we use. We are able 
to accurately asses the disorder strength which places 
the model in the DM phase by considering its finite size 
dependence as p(0) is L independent in the DM phase, 
as shown in Appendix [BI It is possible that the method 
of extrapolation can under- or over- estimate the critical 
disorder strength, and we discuss below in the next sub¬ 
sections the implications of this in determining the criti¬ 
cal exponents. (It turns out that the precise quantitative 
determination of Wc is (is not) crucial for determining 
the critical exponent ly (z) accurately.) 


1. Determining z 

We begin by discussing the procedure we use to esti¬ 
mate the dynamical exponent 2 (see Fig. [3]). After deter¬ 
mining Wc by extrapolation we perform a power law fit to 
the energy dependence of p{E) a,tW = Wc [see Fig.[3][a)] 
and then use the scaling formula (see Appendix A) 

p{E) ~ (10) 

to get z, here we hnd z = 1.46±0.05 for both BD and GD 
(we only show the BD results in Fig. 3 (a)). Error bars 
are estimated by considering the value of 2 due to the 
uncertainty of Wc, which is actually quite small. Even 
arbitrarily assuming a value Wc that exceeds the error 
bars, for example Wc = 2.65 (where the finite size effects 
are still large), yields a value of z close to 1.55 for BD, 






9 


and thus an inaccurate estimate of Wc does not produce 
a large deviation in the extracted value of 2 . We then 
check for consistency in the estimated 2 by computing the 
specific heat, which integrates over the entire bandwidth, 
as shown in Fig. [3jb). Here we use the scaling of the 
specific heat at the QCP 

C'v(T)~T‘^/^ (11) 

and the result from Cy is z = 1.48 ± 0.05 in good agree¬ 
ment with the DOS result shown in Fig. [SKa). We men¬ 
tion that very similar results are obtained for GD which 
are not shown here. The fact that the same value of 2 ; 
is obtained numerically from independent considerations 
of the DOS and specific heat gives high confidence in the 
numerical accuracy of our estimated dynamical exponent. 

2. Determining v 

We now turn to estimating the correlation length expo¬ 
nent V (see Fig.|T|). Here, we find that obtaining a precise 
value of V using the KPM method alone is numerically 
challenging as the value is sensitive to both the accuracy 
of Wc and the range over which the numerical results 
are fitted to extract the power law dependence. These 
features in the KPM data for the zero energy DOS has 
been previously pointed out in Ref. [13 and therefore we 
analyze the data accordingly. We first determine v based 
on the scaling 

p(0,T) - (12) 

for sufficiently large L (i.e. where the data is L inde¬ 
pendent) using the values of Wc and z we have already 
obtained. We do not fit the data sufficiently close to 
the QCP, as there are significant finite size effects as 
shown in Fig. |4] (a) and (b). It is well-known that the 
numerical data close to the QCP suffer from severe finite 
size effects since the correlation length exceeds the sys¬ 
tem size around the QCP. As depicted in Fig. U (a) and 
(b) using the thick dashed line, fitting a range of S for 
L = 130 starting where p{0,L) saturates in L, yields a 
value r = 1.48 ±0.3 for BD and 1.36 ±0.3 for GD, where 
the error bars are obtained by considering the inaccuracy 
of Wc, giving consistent results with varying the range of 
the fit. Fitting to a lar ger range of <5 away from the QCP 
(as done in Refs, [l^ [^) as shown in Fig. |4] (a) and 
(b) (the thin dashed line) yields much lower estimates 
1 / = 1.02 ± 0.2 for BD and 1.15 ± 0.2 for GD. We remark 
that fitting the typical DOS (i.e. the geometric averaged 
local DOS) data versus S in Ref. [l^ suffers from these 
same effects and the critical exponents governing the av¬ 
erage and typical DOS at the QCP remain distinct. 

We also consider the effect of an underestimate of Wc 
on the value of v. Again taking the box distribution 
with Wc = 2.65 for the critical disorder yields a value of 
r = 0.92, which points to a large deviation of r with re¬ 
gards to the accuracy of Wc- Thus, the numerical KPM 


technique, while being excellent for estimating the dy¬ 
namical exponent z, is lacking in accuracy in estimating 
the correlation exponent r, a situation quite common in 
numerics on QPTs where the numerical extraction of r 
is often much more challenging as it is strongly affected 
by finite size and fluctuation effects. 

We now come to finite size scaling and data collapse of 
p{0,L). Here we focus on the general scaling form 

p{0,L) = (13) 

in the vicinity of the QCP. Using the values we have 
extracted for Wc and z, we perform data collapse for 
W > Wc as shown in Fig. H] (c) and (d). We denote the 
value of the correlation length exponent extracted from 
finite size scaling as vl- We find excellent scaling collapse 
of the KPM data for equal to 1.46 ± 0.25 for BD and 
1.42 ± 0.3 for GD. Consistent with the power law fit of 
Eq. (HU, we again find that performing the collapse over 
a larger range away from the QCP or accounting for an 
underestimate of Wc yields a value of closer to 1 (than 
to 1.5). 

The fact that both box and Gaussian distributions give 
consistent results for z/(~ 1.5 or 1 depending on whether 
the scaling fit is being carried out close or away from 
the QCP) within the error bars provides evidence for 
the universality of the critical exponent 1 / (i.e. that r 
is independent of the disorder distribution). However, 
while it is possible to estimate r from the KPM data, it 
seems to suffer large fluctuations from small systematic 
errors (e.g. the precise value of Wc). Thus we make 
the conservative estimate for v to lie within the range 
0.9 — 1.5. Much more computationally demanding work 
on much larger systems would be necessary to obtain an 
accurate estimate of the correlation length exponent v 
in this problem. By contrast, our numerically estimated 
dynamical exponent z is reliable and stable. 

3. Scaling in the guantum critical fan 

To establish the universal nature of the quantum criti¬ 
cal fan we compare the numerically determined crossover 
functions TZ and C with their approximate analytical 
forms (derived in Appendix I3) in Fig. [5l when \E\ > 
L-i§-vz rp ^ Focusing on our numeri¬ 

cal data that is restricted to lie in the non-Fermi liquid 
quantum critical fan (i.e. <g; E/t,T/t <C A/t) we 

plot the numerical results for the DOS versus ES~''^ and 
the specific heat versus T5~''^ using the values of critical 
exponents z = 1.46 and ~ 1-46 for the box disorder as 
well as the cross over functions TZ and C determined ana¬ 
lytically on either side of the QCP in Fig. [5l We find that 
all of our data collapse onto a single curve manifesting 
excellent critical scaling. Even though there is a signifi- 
eant difference between the numerically obtained value of 
zz(Ri 1.4) and the approximate one loop result {v = 1), 
the crossover functions obtained from the two methods 
appear to show considerable agreement. We stress that 
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this is not a fit, and the only parameter that is adjusted 
is the non-universal coefficient of the analytical scaling 
functions. Finally, the consistency between our numeri¬ 
cal calculations and the scaling hypothesis (contained in 
Eqs.IH) has led us to conclude that the QCP is indeed 
non-Gaussian and the disorder-tuned SM-DM QPT in 
Dirac systems cannot be described by a Gaussian theory. 

Motivated by the universality of the non-Gaussian 
QCP, in the following section we construct the effective 
field theory of the long wavelength fluctuations, while ac¬ 
counting for the strong coupling between the underlying 
order parameter and the itinerant fermions. 


IV. ORDER PARAMETER FIELD THEORY 

Within the replica field theory formulation of disor¬ 
dered systems, the sample to sample fluctuations are cap¬ 
tured in terms of a disorder induced effective statistical 
interaction between the fermions although the starting 
bare Hamiltonian is noninteracting [cf. Eqs. ([T]) and 
O]. The collective modes of the density fluctuations 
are captured by a bosonic matrix field Q ^ i'i’ai’l) 
Tr{Q) oc p(0) acts as the order parameter. The devia¬ 
tion of Tr{Q) from its vacuum expectation value consti¬ 
tutes the amplitude or the DOS fluctuations, while the 
transverse part of the matrix field captures the slow diffu- 
sons or Goldstone modes. At the SM-DM QGP, both the 
amplitude Tr{Q) and the transverse diffuson modes are 
gapless and remain strongly coupled with the underlying 
fermions. As we show below, the average DOS controls 
the residue of the diffusion pole, and it therefore vanishes 
as VF —>■ W^. Similarly, the quasiparticle residue of the 
Dirac excitations vanishes as W —>■ W~. Consequently, 
inside the critical fan the entire notion of weakly coupled 
quasiparticle excitations becomes invalid and an emer¬ 
gent strongly coupled non-Fermi liquid state is realized. 
The effective Lagrangian for the semimetal-metal QCP 
is therefore described by 

V'o] = ^Dli’a] + CB[Qab] + iMli’bQab, (14) 

where Co describes the Dirac fermions, £5 is a Landau- 
Ginzburg functional of Q up to quartic order, and A is 
the coupling between the fermions and the bosons. This 
effective theory can be analyzed within a d = 4 — e ex¬ 
pansion around the upper critical dimension of 4, which 
leads to a critical value of A ~ 0(e). It is this strong 
coupling between the collective modes and the fermions 
which unifies the problem at hand with the underlying 
conceptual theme of the itinerant QCP in strange metals. 

For deriving an order parameter theory, which can cap¬ 
ture the low energy properties of both the semimetal 
and the diffusive metal, we consider the product of 
the advanced and retarded propagators G^{E -|- w /2 -|- 
ir], x)G^{E — u;/2 — irj, x). This can capture the diffuson 
or Goldstone modes inside the metallic phase. For this 
correlation function, we can write a generating functional 


e~^, where the effective action is 

S = y -I-T 3 - id]T. (15) 

Here 'I'l = —^( 4 ), and T 3 = diag(l,—1) is a Pauli 
matrix that operates on the retarded and advanced la¬ 
bels. All physical quantities such as the average DOS and 
the density-density correlation functions can be obtained 
from the generating functional after taking derivatives 
with respect to the appropriate source terms. For sim¬ 
plicity we again consider only a single Dirac cone and the 
case of a random scalar potential. After averaging over 
disorder with the replica method we arrive at 

S' = y d^'x 'Si[E - (^ + '^3 + ivajdj]'lia 

-J I ( 16 ) 

When uj = 7 ] = 0, the retarded and the advanced sectors 
appear in a symmetric manner, implying there is a flavor 
symmetry between these two sectors. Therefore, 77 acts as 
an infinitesimal external field for the bilinear and 

for A > Ac (p( 0 ) ^ 0 ) this symmetry is spontaneously 
broken. We perform the following Hubbard-Stratonovich 
transformation 



where Q constitutes a matrix order parameter, and 
{Qab) ~ (4'a'I'j)). By computing the susceptibilities for 
different ordering channels, we find that the most favor¬ 
able choice for Q inside the metallic phase is Q = iQoTs, 
for which r] acts as the external field. For 2 < d < 4 the 
integral in the gap equation for Qq 



The order parameter’s expectation value Qo acts as the 
inverse life-time of the Dirac quasiparticles. Since in 
the mean-field calculation we have not accounted for the 
correction to the real part of the self-energy, which leads 
to a modified dispersion relation at the QCP, the mean- 
field solution is inadequate for capturing the scaling form 
Qo ~ (5^^. We expect that accounting for strong order 
parameter fluctuations will finally lead to the correct re¬ 
sult. 

After the Hubbard-Stratonovich decoupling, if we in¬ 
tegrate out the fermion fields, the fermion bubble con¬ 
tributes to the Qab — Qab correlation function. Since the 
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replica indices for such a correlation function are fixed 
from the outset, the fermion bubble is not proportional to 
Nr and it survives the Nr ^ 0 limit. We can further de¬ 
compose the matrix form of Qab = Qab,^i,sT^ ® Fg, where 
Fg are the sixteen 4x4 matrices operating on the origi¬ 
nal spinor index and are 2 x 2 matrices operating on 
the R/A indices (we have not yet considered the Cooper 
channel). By computing the zero momentum part of the 
fermion bubble, we can show that only the conventional 
and chiral density channels corresponding to Fq = / 4 x 4 
and F 5 = 75 channels are most attractive and degener¬ 
ate. Therefore, only these two density channels can have 
gapless diffuson modes. Since the order parameter or the 
amplitude of Q vanishes at the QCP, we can not safely in¬ 
tegrate out the gapless fermion modes. For this reason, it 
is necessary to (phenomenologically) introduce a fermion- 
diffuson coupling (of the Yukawa form) to capture the 
interplay of gapless order parameter and fermionic exci¬ 
tations. For simplicity if we consider E = uj = 0, the 
effective order parameter field theory takes the following 
form 


which implies rjQ = 2z -\- 2 — d (putting in our numerical 
estimate of z yields rjQ za 1.9). We have also introduced 
the diffusion constant (D) 

D - ~ ^-Az+2-d)^ ^26) 

Upon approaching the QCP in c? = 3 from the metallic 
side, the diffusion constant diverges as - 1 ) M . But 
most importantly, the residue of the diffusion pole van¬ 
ishes. Therefore, the diffusion pole loses its meaning as 
a well defined excitation, when the QCP is approached 
from the metallic side. In addition, such an order pa¬ 
rameter theory also gives rise to an anomalous dimension 
for the fermion field ^ 0(4 — d), which implies that 
the quasiparticle residue of the Dirac fermion vanishes as 
when the QCP is approached from the SM side. 

The vanishing of the quasiparticle residue of the Dirac 
fermions can also be addressed within the 2 -|- e expansion 
scheme. After using the one loop RG procedure as in 
Ref. Q the retarded propagator in the SM phase acquires 
the form 


£['Sa,Qab] 

£B[Qab\ 

£DB\f^ah^ 


+ LslQab] + FDsiQab, 'I'a],(20) 

iv I d‘^x 'I'IjQ: • Vika, (21) 



iTr(VQtvQ) + ^Tr(Qtg) 


+wiTr[(Q+Q)2]+M2[Tr(QtQ)]2 


( 22 ) 



d’^x 


(23) 


where the summation over repeated replica indices is im¬ 
plied. We are only considering Qab in the regular and 
the axial density channels. When an external frequency 
is considered, it couples to Tr[Qr 3 ] in the regular density 
channel as an external field. 

At the critical point r = 0, both regular and axial 
density fluctuations become gapless. If we scale x ^ xe\ 
4- Q find A 

Therefore, d = 4 serves as the upper 
critical dimension, and the interaction effects can be ad¬ 
dressed through a d = 4 — e expansion, which leads to a 
critical point at Ac = 0 (e), uj = 0 (e). 

In the DM phase the density-density correlation func¬ 
tion (involving both retarded and advanced sectors) dis¬ 
plays a diffusion pole corresponding to the gapless Gold- 
stone modes or the diffusons, which behaves as 


nd(w,k) 


Q\ 


iQow - S-'^Q 


ILJ 


- Dk‘^ 


(24) 


where tjq is the anomalous dimension of the gapless col¬ 
lective mode. If the order parameter inside the DM phase 
varies as Qq ^ 6^ , the hyperscaling relation leads to 

13 = {d - 2 + tiq)^ = ziy, (25) 


Gr{E = 0) 


[r;Fo -f vSG-£''k ■ r] ■ 


(27) 


Therefore the quasiparticle residue and the effective 
Fermi velocity are respectively given by and 

vS'^G-i)^ At the one loop order = (z — I) = e/2 and 
iz = 1/e. Consequently, both the quasiparticle residue 
and the Fermi velocity behave as ^ This is an arti¬ 
fact of the one loop analysis and at higher loop orders 
(beginning at e^) they will follow different power laws. 
Hence, we conclude that the quasiparticle residue and 
the effective Fermi velocity of the Dirac fermion vanish, 
as we approach the QCP from the semimetal side. 

Consequently, the residue of the two different types of 
quasiparticles vanish while approaching the QCP from 
either side. Thus, the critical region of this system can 
not have any simple quasiparticle description, qualita¬ 
tively similar to what is envisaged for QC regions in cor¬ 
related materials. Thus, the orange region in our nu¬ 
merically obtained quantum phase diagram of Fig. [1] is 
a finite-temperature ‘non-Fermi liquid’ quantum critical 
crossover regime arising from the critical fluctuations of 
the non-Gaussian QCP underlying the SM-DM QPT in 
the three dimensional Dirac materials. 

To summarize this section, the finite expectation value 
of the DOS [as shown numerically in Fig. [2l[a)] on the 
metallic side and its critical fluctuations in the vicinity 
of the QCP can be described by an order parameter in the 
form of the lagrangian Lb for the matrix field Q. While 
approaching the QCP from the SM phase, the quasiparti¬ 
cle residue of the underlying Dirac fermions (as described 
by Lb) vanishes continuously as (5 —>■ 0. Thus, the gapless 
fermionic and bosonic degrees of freedom must be treated 
on an equal footing and this has led us to construct the 
effective action in Eqs. (nnn-dMi), which governs all the 
long wavelength degrees of freedom. This effective the¬ 
ory retains both longitudinal and transverse fluctuations 
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of the Q matrix, as required by the vanishing of the av¬ 
erage DOS at the QCP. This should be contrasted with 
the familiar nonlinear sigma model description of the An¬ 
derson localization at stronger disorder [I^, where the 
average DOS remains noncritical across the transition. 
The nonlinear sigma model is obtained by integrating 
out the gapped ballistic fermions and the gapped longi¬ 
tudinal fluctuations of the Q matrix inside the metallic 
phase (below the energy scale of ^ Qo)- 

In our earlier work dl on the phase diagram of a dis¬ 
ordered DSM, we have noticed a significant difference 
between the critical exponents for the typical DOS at 
the SM-DM QCP and the QCP describing an Anderson 
localization transition (at stronger disorder). We believe 
that the differences between the multifractal properties 
at these two QCPs are inherited from (i) the presence or 
absence of itinerant, ballistic fermions, and (ii) the pres¬ 
ence or absence of gapless longitudinal mode of the Q 
matrix. A detailed analysis of different correlation func¬ 
tions and the multifractal properties associated with the 
field theory in Eqs. dinii-disi) is beyond the scope of this 
work and remains an open problem for the future. The 
Dirac SM and the DM phases preserve time reversal and 
inversion symmetries, and do not display any topologi¬ 
cal transport properties such as anomalous Hall effect. 
For this reason, there is no Chern Simons term in the 
effective field theory. By contrast, a time reversal sym¬ 
metry breaking Weyl semimetal phase and the resulting 
diffusive metal possess an anomalous Hall effect. In the 
effective field theory this arises due to a topological Chern 
Simons term [s^. 

One direct implication of our effective field theory is 
the existence of an upper critical dimension of du = 4 
for the disorder-tuned SM-DM transition in Dirac-Weyl 
systems (in addition to the known lower critical dimen¬ 
sionality of d; =2) . For dimensions d > d„, the correla¬ 
tion length exponent v should be given by the mean-field 
value iz = 1/2 and hyperscaling relations will be vio¬ 
lated. It will be interesting to check this observation in 
future numerical studies on higher dimensional models 
of a Dirac SM. Finally, the existence of an upper critical 
dimension for the SM to DM QCP will further distin¬ 
guish this transition from Anderson localization, which 
does not possess an upper critical dimension. 

We mention here that we do not, however, have an ex¬ 
planation for our numerical finding of the exact z(~ 1.5) 
being equal within the numerical accuracy to the one- 
loop theoretical value (z = d/2). Whether this is a mere 
coincidence or a deeper truth (i.e. somehow all higher- 
loop corrections cancel out) is unknown at this stage. Of 
course, our numerical error in estimating z, while being 
small, is not zero, and thus small corrections in the ex¬ 
act theoretical value of z (well less than 10%) away from 
z = 1.5 cannot be ruled out without more numerical work 
with much larger systems. There are well-known exam¬ 
ples in the literature for numerical and one-loop theoret¬ 
ical z values being very close to each other purely for¬ 
tuitously in completely different contexts of dynamical 


phase transitions 


V. DISCUSSION AND CONCLUSION 

The experimentally observed non-Fermi liquid scaling 
behavior of many correlated metals over a sig¬ 

nificant range of temperatures is usually reconciled with 
the existence of a large quantum critical fan driven by 
the putative QCP, whose intrinsic properties are often 
not well-understood. It is generally believed that the 
generic behavior of the ‘strange metal’ within the finite- 
temperature critical fan regime is non-Fermi liquid like 
because of strong quantum fluctuations arising from the 
QCP (which may often be ‘hidden’ or ‘inaccessible’ ex¬ 
perimentally), but no general proof exists. For describ¬ 
ing such a metallic or itinerant QCP, it is necessary to 
address both the gapless fermionic degrees of freedom 
and the bosonic order parameter fluctuations on an equal 
footing, which is a highly challenging task for analyti¬ 
cal calculations, particularly for strongly correlated sys¬ 
tems, where even the minimal model is intractable with 
or without the QCP. In this work even though we have 
focused on non-interacting massless Dirac fermions we 
have been able to make various connections to the general 
notion of scaling phenomena of itinerant QCPs. First, 
we find the appearance of a broad quantum critical fan 
at finite temperatures. Second, we have succeeded in 
constructing an effective field theory of bosonic fluctua¬ 
tions (captured here by the Q matrix for diffusive modes) 
that are strongly coupled to itinerant electron degrees of 
freedom (i.e. massless Dirac fermions). Third, we have 
shown how the general notion of quasiparticle excitations 
breaks down (via their residue) in the quantum critical 
fan region even for our noninteracting model, explicitly 
demonstrating that a QCP by itself, independent of the 
underlying model being interacting or not, leads to a 
generic non-Fermi liquid behavior. However, as we are 
dealing with non-interacting fermions, it is important to 
contrast the version of itinerant quantum criticality we 
obtain here with the form that is applicable to strongly 
correlated electronic systems. We are dealing with free 
fermions and the interaction is statistical and only due 
to disorder (only elastic scattering). Consequently, the 
effective field theory will not be a dynamical one (with 
inelastic scattering) as we have found. This is in sharp 
contrast to that of the strongly correlated systems where 
the dynamics of the bosonic field is also important and 
must be incorporated in the theory on an equal footing 
with the fermionic field. 

Through a precise numerically exact calculation using 
the KPM technique and an approximate field theoretical 
analysis, we have established that the three-dimensional 
Dirac SM phase undergoes a disorder-tuned QPT into a 
DM state for different types of disorder (potential, axial, 
and mass), which belong to the same universality class 
(and surprisingly, even with the identical critical coupling 
for different disorder types) of an itinerant QCP. We find 
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the critical exponents for box and Gaussian distributions 
agree to within numerical accuracy (with z = 1.46 ±0.05 
and V £ [0.9,1.5]), which points to the universal nature 
of the QCP. For this model of an itinerant QCP, the 
universal scaling functions obtained from numerical and 
analytical methods show a remarkable agreement, and 
thus can serve as a conceptual framework for address¬ 
ing other itinerant critical phenomena in diverse corre¬ 
lated metallic systems. Experimentally, a dilute three- 
dimensional Dirac material such as NaaBi with a small 
value of the Fermi energy Ep A associated with resid¬ 
ual doping, can be a promising material for verifying our 
theoretical results. In a real system, the residual uninten¬ 
tional doping always produces a finite Fermi energy shift¬ 
ing the system from the zero chemical potential Dirac 
point, but as long as this doping-induced Fermi energy is 
lower than the temperature, one expects intrinsic Dirac 
point behavior to prevail, leading to the manifestation 
of the QPT studied in this work. One experiment could 
be the measurement of the specific heat in the quantum 
critical regime which, according to our theory, should 
manifest a non-Fermi liquid like quadratic temperature 
dependence associated with the critical quantum fluctu¬ 
ations. One key aspect of this QCP is its fundamental 
non-Gaussian nature. Our KPM-based estimate of the 
correlation length exponent v, however, has rather large 
error bars, and much more work would be necessary (well 
beyond the scope of the current work) to obtain a pre¬ 
cisely numerically accurate value of ly in the SM-DM QPT 
in three dimensions. 

Finally, we comment on the possible role that ‘Grif¬ 
fiths physics’ (associated with rare fluctuations) might 
play in the SM-DM QPT studied in the current work. 
The physics of rare regions (or the so-called “Griffiths 
physics”) is not captured by the long wavelength field 
theory developed in the present manuscript. Recently it 
has been suggested that rare regions affecting the SM- 
DM QCP can induce a finite DOS at the Dirac point 
for an infinitesimally weak disorder [d^, thus converting 
the SM phase into effectively a DM phase. Compared to 
the box distribution, the unbounded tails of the Gaussian 
disorder distribution are expected to enhance these rare 
fluctuations, but we get the same universal QCP prop¬ 
erties for both box and Gaussian distributions without 
seeing any obvious effects of rare regions. Within our 
numerical calculations we have not found any evidence 
of rare region effects, as the results for both distribu¬ 
tions agree with the field theoretical predictions for the 
existence of the transition (e.g. considering the problem 
in two versus three dimensions). However, the detec¬ 
tion of rare region effects can be quite subtle (due to a 
possibly exponentially small DOS contribution from the 
rare regions) and therefore the possible existence of rare 
regions still cannot be ruled out completely (e.g. the re¬ 
gions are rarer than our numerical system size or the rare 
regions being independent of the QCP itself or the expo¬ 
nentially small DOS introduced by the rare regions are 
simply too small to show up in our KPM simulations). 


which necessitates a detailed focused numerical study on 
its own. Although we cannot rule out the possibility of 
rare regions based on our study, we can assert that the 
disorder-tuned SM-DM QPT is well-captured by our the¬ 
ory and numerics once any contribution from rare fluctu¬ 
ations are subtracted out. It is important in this context 
to emphasize that the real Dirac systems are likely to 
have random impurity-induced ‘puddles’ or spatial den¬ 
sity inhomogeneities (i.e. random spatially nonuniform 
doping) around zero energy (i.e. the Dirac point) any 
way, as is ubiquitous in graphene [1^, leading to a lo¬ 
cally fluctuating chemical potential masking the QCP in 
a way similar to what rare regions would do. Any exper¬ 
imental study of the QPT itself must find a reasonable 
way of subtracting out these additional effects (from rare 
regions and/or puddles) in order to capture the underly¬ 
ing critical behavior. 
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Appendix A: Scaling formulas 

In the vicinity of the QCP, we have two divergent 
scales, namely the spatial correlation length oc 6 ~''' 
and the temporal correlation length ~ oc The 

universal scaling functions are determined by these two 
divergent scales. Quite generally for an observable O we 
have the scaling formula, 

0{L,E,T) ^ (Al) 

where L, E and T respectively denote the system size, 
the energy and the temperature, and Aq is the scaling di¬ 
mension of O. For example, the scaling dimensions Aq 
for the average density of states, the free energy den¬ 
sity, the specific heat and the longitudinal conductivity 
are respectively given by —(d — z), —(d ± z), —d, and 
— (d— 2). The leading scaling behavior is determined by 
max{L“^, A, T, Since there are two fixed points 

respectively with z = I and z = 3/2 (at one loop), the 
crossover scaling functions are quite nontrivial. 

For example, first consider the average density of states 

piL,E) ^ ACf). (A2) 

At zero energy this becomes 

p{L, 0) - ~ (A3) 
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When L ^ S , we are in the quantum critical 

regime and 

(A4) 

reflects the scale invariance. By contrast, inside the 
metallic phase (A > A^) 

p{L,0) ^ (A5) 

On the semimetal side we rather expect 


where Aq = A(/ = 0) is the bare value of the coupling 
constant. On the SM side Aq < A^, the renormalized 
disorder coupling monotonically decreases and satisfies 
the asymptotic form 

A(0 (A14) 

0 

In contrast, for the metallic side Aq > Ac the renormal¬ 
ized disorder coupling diverges at the RG scale 

e'l =e^A^|d|-^ (A15) 


(A6) 

These finite size scaling properties constitute the basis 
for the analysis of our numerical results. In addition, 
we consider a large enough system size, L :s> {^i,v/E). 
In this case, the scaling behavior is determined by the 
interplay between v/E and . The scaling function now 
behaves as 


p{E) - (A7) 

When, E ^ we are in the quantum critical regime 
and obtain the following scale invariant answer 

p{E) - (A8) 

Inside the metallic phase (A > Ac), the leading behavior 
is given by 

p{E) ^ . (A9) 

By contrast, inside the Dirac semimetal phase we obtain 

piE) ~ c{5)E‘^-^ ~ (AlO) 

We can gain considerable analytical insight regarding 
the crossover functions by integrating the one-loop renor¬ 
malization group equations. The one loop RG equations 
for potential and axial disorder are @ 

zil) = 1 -b A, = -(d - 2)A -b 2A2. (All) 

dl 

We simultaneously solve the beta function for the disor¬ 
der coupling A and 

^ = z(0A = (l + A)e, (A12) 

where e = E/{vA) is the dimensionless quantity de¬ 
fined from the energy E. The dimensionless tempera¬ 
ture t = kBT/{vA) also satisfies the same flow equation 
as e. From these solutions we can identify two RG flow 
times associated with two divergent scales. The disorder 
coupling A(/) is given by 


which defines the correlation length and the scaling 
exponent v = l/(d— 2). By contrast, the divergent time 
scale has to be found from the solution of the differential 
equation for energy 


e(0 = e(0) 



A(0-AcyY Ao \ 

A{0)-aJ \Ail)) 


e(0)e^' 




(A16) 


For that we need to find the RG scale I 2 as a function 
of the bare energy e(0) and <5 by imposing the condition 
^ih) ^ 1- This results into the following equation for I 2 






1 -b 




iz-l)u 


(A17) 

When A = Ac, Zi —>■ 00, and the infrared cutoff is solely 
determined by I 2 where 




Inside the semimetal phase A —>■ 0, and again the infrared 
cutoff is determined by I 2 with ^ . Since p 

has the scaling dimension —(d—z), we find p ^ 
inside the critical fan, and p ~ c{S)\E\^‘^~^'> inside the 
semimetal phase, with the c(S) quoted in the previous 
paragraph. 

The explicit solution for in d = 3 inside the quan¬ 
tum critical fan is obtained by finding the roots of a cubic 
equation. When Aq > Ac, we have the following cubic 
equation 




-e ^ — 


l + S 


= 0 , 


£ 2 ( 0 )'' £ 2 ( 0 ) 
which has the only one real root given by 

1 . . -1 I (1-b (5)3V3e(0) 


(A18) 


= 2, 


3e2(0) 


sinh 


■ sinh 


2 ( 53/2 


36A 


sinh 


1 


sinh ^ a; I , for (5 << I, 


(AI9) 


A(0 


I-b 



g(d-2)/ ’ 


(A13) with the anticipated dimensionless variable x = 
——M-. When X >> 1 we have the quantum critical 
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behavior, which eventually gives away to the metallic be¬ 
havior for X < 1 when the infrared scale is determined by 
or the correlation length. When Aq < Ac, the cubic 
equation becomes 


g3Z2 _ 


1-5 


:e — 


£ 2 ( 0 ) £ 2 ( 0 ) 


= 0 . 


(A20) 


When 1 


When \E\ » , these analytically obtained 

approximate crossover scaling functions have been com¬ 
pared with the numerically determined exact scaling 
function TZ in section uni In this comparison, we have 
only adjusted an overall numerical prefactor for the an¬ 
alytical formula. A similar approximate expression for C 
is also found by substituting the explicit expressions of p 
(in terms of g{x)) in Eq. [51 


= 2 , 


3£2(0) 


cosh 


i cosh 

O 


-1 


cosh (- cosh ^ X 
X V 3 


(1 - 5)3y3£(0)' 

I 2^3/2 

for,5«l. (A21) 


For ^ there are three real roots, and only 

one of them is positive. This root is given by 


e’'^ = 2. 


3£2(0) 


36A 


cos 


1 _1 j (1 - 5)3V3£(0) 

3 I 2^3/2 


cos ( - cos X , for (5 << 1, 


(A22) 


where the inverse trigonometric function is restricted 
to the first quadrant. For small 5, 5e^'^ is again a 
scaling function of £(0)(5“3/2_ Por energies satisfying 
£(0)(5“3/2 >> 1 or X >> 1 inside the quantum critical 
fan cosh and sinh terms determine the z = 3/2 quan¬ 
tum critical behavior, while inside the SM phase the cos 
term leads to the z = 1 critical behavior. Therefore, the 
crossover from the 2 = 3/2 to 2 = 1 scaling is appropri¬ 
ately captured by how the cosh term transforms into a 
cos term. 

In the thermodynamic limit the total number of states 
per unit volume can be estimated to be 

® ~ = ^r‘^f-‘^ix). (A23) 

1 j^ 

Since the average DOS can be obtained by differentiating 
N{e)L~‘^ with respect to £, we can write 

. ^ A 253/2 ^ ^ ^ ^ ^ 

P(e) ~ 2V3 u ^ ^ 

where the crossover scaling function g{x) is determined 

X coth (i sinh”^ x) 

3-\/l -|- x2 

(A25) 

X tanh (i cosh”^ x) 
3-\/x2 — 1 

(A26) 

X tan (i cos“3 x) 

3-\/l — x2 

(A27) 


by 

aix) 


9{x) = 


9{x) = 


cosh® (i cosh ®x) 

for A < Ac, X > 1 
.2 


X 


=3 ( 1 , 


3 ■ 3 ) 

for A < Ac, X < 1. 


1 
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FIG. 6: (a) (Color Online) Finite size scaling of p(0, L, Nc) 
on a log-log scale, for Nc = 1028 in the clean limit {W = 0) 
yields a power law decay in system size consistent 

with IfL^. (b) Scaling of p{0, L, Nc) on a log-log scale with 
Nc for L — 60 yielding a linear in Nc scaling. Here Nc is 
the number of terms kept in the Chebyshev expansion in the 
KPM numerics. 




Appendix B: Numerical details 

In this appendix we discuss the scaling and conver¬ 
gence of the zero energy DOS as a function of various 
KPM parameters. We begin with the clean limit W = 0 
and compute the zero energy density of states as a func¬ 
tion of system size ranging from A = 60 to 140 in steps 
of 10. Generally in the thermodynamic limit we expect 
p{0,L) ~ 1/^2. However, this behavior cannot be ob- 
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FIG. 7: (Color Online) DOS for L = 60 and IF = 2.01 < Wc 
(a) dependence of p(0, L, N^) converges in for > 
4096. (b) Low energy dependence of p(E, L, Nc) for various 
Nc showing that the deviations in the DOS for larger Nc in 
the SM phase is a result of the KPM beginning to resolve 
individual eigenstates. The low energy peak is a result of the 
8 Dirac states that are there in the clean limit that have now 
been broadened by disorder. 


gree of accuracy as shown in Fig. [6Kb). These results are 
obtained by considering the asymptotic behavior of the 
Jacobi theta function 

OO 

= 03 ( 0 , e-^^) = t?(0, ix). (B2) 



9000 12000 15000 18000 

N. 



tained within the KPM. Rather within the KPM we 
find p(0,L) ~ consistent with a 1/L^ scaling 

[see Fig. EKa)]. This can be understood by considering 
the density of states with a Gaussian broadening fac¬ 
tor a = tt/Nc, where Nc is the number of terms kept 
in the Chebyshev expansion (as described in section Hill 
A), as this is an excellent approximation to the effect of 
the Jackson kernel [j^ . Since there is only an intensive 
number (8 corresponding to each cone) of states at zero 
energy the KPM yields 


p{0,L,Nc) = 


exp 


- + nl) 


rix ,ny ,nz 

1 Nc 

crL^ ttL^ 


L^V2Tra‘^ 


(Bl) 


where we have used ki = 2TmilL. We find p(0, L, Nc) ^ 
reproducing the linear in Nc scaling to a good de¬ 


FIG. 8: (Golor Online) DOS for L = 60 and IF = 3.0t > Wc, 
(a) Nc dependence of p(0, L, Nc) converges in Nc for Nc > 
4096. (b) Low energy dependence of p{E, L, Nc) for various 
Nc showing that there is a very weak dependence in the DM 
phase. 

Moving away from the clean limit by introducing dis¬ 
order provides a natural broadening of the energy eigen¬ 
values and thus changes the scaling with Nc- Focusing 
on the SM phase, we find that p{0,L,Nc) Xjl? is 
essentially satisfied for weak disorder strengths. Even 
though we find a large deviation in the zero energy DOS 
as a function of Nc, it does converge for increasing Nc as 
shown in Figs.|7|(a) and (b). The deviations in p(0, L, A^c) 
are due to the KPM beginning to resolve the energy 
eigenstates. 

Inside the DM phase, p(0, L, A^c) has a very weak de¬ 
pendence on Nc [as shown in Figs. |5| (a) and (b)] and for 
sufficiently large system sizes, converges in L [see Figs.lSj 
(a) and (b)[. We conclude that in each phase choosing 



















Nc = 4096 gives a converged estimate of p(0), and there¬ 
fore we use this value of to determine the critical 
exponents. 
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(b) 0.01 



0.001 


0.0001 



FIG. 9: (Color Online) DOS at zero energy p(0) as a function 
of of the inverse system size 1/L on a log-log scale in both 
the SM and DM phases for potential disorder up to L = 
110. Results for the box distribution (a) and the Gaussian 
distribution (b). In the SM phase p(0) goes to zero as the 
system size is increased, while it saturates in the DM phase, 
and at the QCP p(0) is still decreasing for increasing L. In 
the DM phase the results are saturating in 1/L. 
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